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Abstract 

We present a lattice NRQCD study of the B meson decay constant in the 
quenched approximation with emphasis given to the scaling behavior. The 
NRQCD action and the heavy-hght axial current we use include all terms 
of order 1/M and the perturbative ©(a^a) and 0{as/M) corrections. Using 
simulations at three value of couplings /?=5.7, 5.9 and 6.1 on lattices of size 
12^ X 32, 16^ X 48 and 24^ x 64, we find no significant o dependence in fs if the 
0{asa) correction is included in the axial current. We obtain fs = 167(7) (15) 
MeV, /b, = 191(4)(17)(;^) MeV and feJfB = 1.15(3)(l)(t|;), with the first 
error being statistical, the second systematic, and the third due to uncertainty 
of strange quark mass, while quenching errors being not included. 

PACS number(s): 12.38.Gc, 12.39.Hg, 13.20.He, 14.40.Nd 
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I. INTRODUCTION 



Lattice QCD provides a promising approach for a first-principles calculation of the 
hadronic matrix elements of B meson relevant for a precision determination of the Cabibbo- 
Kobayashi-Maskawa (CKM) matrix elements. The most important matrix elements is the B 
meson leptonic decay constant fs- it is needed to determine Vtd] from the lattice technical 
point of view it is the simplest B meson matrix element calculable in lattice QCD with 
which one can study systematic errors associated with a lattice treatment of heavy quark. 

The need for a careful examination of systematic errors stems from the fact that their 
magnitude for naive quark actions such as the Wilson action is of 0{aM) with M the heavy 
quark mass. Hence errors of this origin can exceed 100% for a typical lattice spacing of 
~ 2GeV used in present simulations. To overcome this problem, recent lattice studies 
of J'b m employ non-relativistic effective theory of QCD (NRQCD) or a non-relativistic 
interpretation of the relativistic lattice quark action for heavy quarks 

NRQCD is an effective theory formulated as an expansion in D/M where D is the spatial 
covariant derivative which is of O(Aqcd) in the heavy-light system. For NRQCD one has 
to choose the coefficients of the expansion by imposing a matching condition with the full 
theory. This can be made using perturbation theory. In practice one has to truncate both 
non-relativistic expansion and perturbative expansion at some order so that the systematic 
error in NRQCD calculations is organized as a double expansion in Aqcd/^ and the strong 
coupling constant a^. 

An additional source of systematic errors is the discretization error proportional to some 
power of oAqcd- Since NRQCD is valid only when aM > 0(1), the continuum limit 
a — > can not be taken. Therefore, removing discretization errors is more important in this 
formalism than in the usual relativistic formulations for which continuum extrapolations can 
in principle be made. For this reason, in many lattice NRQCD calculations, the correction 
terms to remove oAqcd and even (oAqcd)^ errors were introduced. 

Until recently the matching coefficients for the action [^|-|^ and the current opera- 
tors [0 have been available only at one-loop level without operator mixing. This means 
that 0{asAQCY)/M) and O^asaAQcr,) errors had been left unremoved. Recently, Shigemitsu 
and Morningstar carried out a one- loop calculation necessary for an 0(q;sAqcd/^) and 
O(asaAQCD) improvement of the heavy-light axial vector current [^0. The first simulation 
including this improvement was performed by Ali Khan et al. [jlO|,|lT|, in which they pointed 
out that the 0{asAQCB/M) and O^asaAqcT)) terms significantly affect the values of fs- 

The study of Ali Khan et al. pO|JTl[] was made at a single lattice spacing corresponding 
to the inverse gauge coupling (3 = G/g"^ = 6.0, and hence left open the important question of 
the lattice spacing dependence of fs obtained with lattice NRQCD (in Refs. |]T2| Hein has 
calculated fs^ at j3=5.7 and discuss the scaling behavior by combining the result at /?=6.0 
of Ref. |TT[). In this article we report on results of our systematic study concerning this 
question. Our simulations are carried out with the plaquette action for gluons at /? = 5.7, 5.9 
and 6.1 corresponding to the range of lattice spacing a ~ 0.18 — 0.09fm. For light quark 
we employ the 0(a)-improved Wilson (clover) action |T^ with the tadpole improved one- 
loop value for the clover coefficient ||l^ , [l5[| . We investigate in detail the effect of one- loop 
improvement of the heavy-light axial vector current as a function of the lattice spacing. 
Our main results are obtained to 0(1/M), but the question of higher order corrections are 
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examined by comparing the results for the NRQCD action complete to 0(1/M^). 

This paper is organized as follows. In Sec. |l| we summarize the NRQCD action we use. 
In Sec. JTT| improvement of the axial vector current is discussed, and our one-loop mixing 
coefficients are presented. Details of the simulations and our methods for extraction of the 
decay constant are given in Sec. |rV| together with numerical results. We discuss the effect 
of improvement in the static limit in Sec. 0. Our results for fs are presented in Sec. |V1 
where a comparison is also made with those obtained with the relativistic formalism. In 
Sec. [VI I| the hyperfine splitting of the B meson and the Bg — B mass difference are given. 
Our conclusions are summarized in Sec. IVTTT. 



II. LATTICE NRQCD ACTION 
A. Form of action 

Let us denote by Q{t, x) the two-component heavy quark field. This field evolves in the 
time direction according to the action, 

S = J2 Q\t, x) [Q{t, x) - KtQ{t -l,x)]. (1) 

t,x 

where the operator Kt specifies the evolution; our choice is 

aiJoA" A aSH\ ^ ^ a5H\ aHo 



2n Jt \ 2 \ 2 J^_^ V 2n 



t-i 



Here subscripts represent the time slice at which Hamiltonian operators such as {l — aHQ/2n) 
act, and an integer n is introduced to avoid instability appearing in the evolution equation 
due to unphysical momentum modes f^. We note that the ordering of terms in Eq. (|]) is 



different from the one employed in |11|: the factor {l—a6H/2) is placed inside of {l—aHQ/2n) 
in our choice. 

The leading order Hamiltonian Hq is given by 

For the correction term 6H, we consider two choices corresponding to the non-relativistic 
expansion to order 1/M {SHi) or to order [SHu), given by 



-C3 
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mi 



cr ■ fA^^) xE-Ex A^^^l 
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We refer to the two choices as NRQCD-I and NRQCD-II. We work with both Hamiltonians 
in parallel and compare their results in order to examine effects of truncation in the 1/M 
expansion. Various covariant differential operators in the Hamiltonian are defined in terms 
of the forward and backward derivatives A^^^ and A^"-* in the fx-th direction as A^''^-' = 

(A(+) + A(-))/2, A(f) ^ A(+)A(-), A(2) = ELi A^, and ^ Eii(Af )2. The field 
strength operators B and E are constructed with the clover-leaf definition as in Ref . 0| . The 
bare heavy quark mass is denoted as Mq, and Cj's are parameters to describe the strength 
of each term. 

The relativistic four-component field iph is related to the effective field Q through the 
Foldy-Wouthuysen-Tani (FWT) transformation, 

^ljh{t,x)=RQ{t,x), (6) 

Here the transformation operator R is given by 

Rji = 1 — di h do 7T 

" ^ 2Mo 8Mi 

with S-' = diaglo"-', a-'}, and Ri {Ru) is to be used in conjunction with 6Hi {6Hii) to achieve 
the desired accuracy in the 1/M expansion. 

The coefficients q and di should be determined by matching the action to the continuum 
relativistic QCD action by either resorting to perturbation theory or estimating it non- 
perturbatively so as to reproduce the same theory in each order of the 1/M expansion. So 
far even perturbative results are not available for these coefficients. We adopt the tree- 
level value Cj = 1 and = 1 in our work, applying, however, the mean-field improvement 
to all the link variables in the action and the FWT transformation with the replacement 
U^/uq, where we take Mq = (Tr?7piaq/3)^/^ |TH. 



III. IMPROVEMENT OF THE CURRENT 

To calculate the decay constant J'b, the heavy-light axial vector current in lattice NRQCD 
has to be matched to that in continuum QCD. For the overall renormalization factor Za this 
was first performed by Davies and Thacker by perturbation theory to one-loop order. 
An important recent development made by Shigemitsu and Morningstar is that the 
matching has been extended to O(asaAQCD) and ©(a^AqcD/-^)- Since our choice of the 
action is slightly different from that used by the authors of Refs. [0,^, we have repeated the 
one-loop calculation for our NRQCD action. 

Consider the axial vector current A^^^^^ in the continuum. We demand that on-shell 
S matrix elements of the lattice axial current reproduce that of the continuum current up 
to 0{p) with p the spatial momentum of the heavy or light quark. At one- loop level the 
relation takes the form 
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4 -\^ ^n, -U n, n^^^ ^^^^ -U n^^) 7(2) ,'n^ 
^4cont -[-■- + ^sPa \ -'latt + "sPa J\a.tt + "sPa -'latt' 

where the heavy-hght lattice operators of dimension 3 and 4 are defined by 

JiSl = v^^rv^,, (10) 

J^l = ^^^r7 ■ A(±)^,, (11) 

i'l = ^^/7 ■ A^r^,, (12) 

with r = 7574 for the temporal axial vector current, and tpi and tph denoting the light and 

heavy quark fields, respectively. We calculate the coefficients p^'' for NRQCD-I for heavy 



quark and the 0(a)-improved clover action ||T^ for the light quark. The use of clover action 
for the light quark is necessary to achieve the accuracy of ©(a^a) in matching the current. 
For renormalization of the continuum current we adopt the MS scheme using dimensional 
regularization with fully anti-commuting 75. We also apply the tadpole improvement pro- 



cedure |T6| with the average plaquette for all link variables in the covariant derivative of 
the operators, ([TT| ) and (0), and with the critical hopping parameter for the wave function 
renormalization of the light quark fields consistently in both non-perturbative and pertur- 
bative calculations. 

Numerical results for the coefficients p^^ are listed in Table |, and plotted in Figure |l] as a 
function of I/oMq. For p^^ the difference p^"* — (l/vr) In(aMo) is shown since the coefficient 
of the leading operator 1 + OsPa'*; being the usual renormalization factor Za, contains a 
logarithmic term (l/vr) In(aMo). The other coefficients p^"* and p^"* are divided by 2aMo. 
The filled symbols represent the values explicitly obtained with the static action We 



have confirmed that the infinite mass limit of - (l/vr) In(aMo) agrees with the static 
results of Borrelli and Pittori |T^ and of Golden and Hill |T9| . 

We observe that p^''/2aMo vanishes in the limit aM 00, which tells us that the con- 
tribution of OsPa^ Ji^tt is of 0{asAQCD/M). This is expected since Ji^^^ involves a derivative 

(2) (2) 

of the heavy quark field. On the other hand, OsPa Jiatt does not contain such a derivative, 
and pi^ /2aMQ remains finite in the static limit as seen in Figure Thus its contribution 
contains terms of ©(a^aAQCD)- This term is an analogue of the current improvement term 
of O^asO) for the light quark discussed in Ref. |T5| . 



We add a remark that we have repeated the one-loop calculation for the action employed 
in Ref. 0, and numerically confirmed their results to a 3 digit accuracy. 



IV. DETAILS OF THE SIMULATION 
A. Run Parameters 

We list our simulation parameters in Table ||. Our simulations are carried out for three 
values of the coupling /3=5.7, 5.9 and 6.1 using the standard plaquette action for gluons. The 
corresponding values of the lattice spacing a is about 0.18, 0.13 and 0.09 fm, respectively, 
if determined from the string tension. We choose our spatial lattice size to be larger than 2 
fm. 
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For heavy quark we take five values of tlie bare mass aMo for each /3 to cover a range 
of the physical heavy quark mass M between 2 GeV and 16 GeV. This wide range of heavy 
quark mass enables us to examine explicitly the 1/M dependence of fs- The parameter n 
is chosen so as to satisfy the stability condition n > S/oMq. 

For light quark we use the 0(a)-improved Wilson action with the clover coeffi- 
cient Csw = (l/^o)[l + 0.199a;y(l/a)], which includes the 0{as) correction calculated in 
Refs. [0,^. Four values of the light quark hopping parameter k are employed for extrap- 
olation to the chiral limit (see Table |I| for numerical values). 

In the quenched approximation the value of the strange quark mass differs depending 
on whether mx or is used as input. The value of determined with is higher than 
that with rriK, and the discrepancy does not become smaller for smaller lattice spacings. 
We choose to calculate for both m^, and take their difference as a systematic error. The 
hopping parameters (k^i from and from m^) are also given in Table ||. 

The physical scale of lattice spacing a is fixed using the string tension a = 427 MeV. Re- 



cent data of the string tension for the standard plaquette action are summarized in Ref. 
We adopt their parameterization to obtain the values of 1/a at our j3. 



B. Fitting procedure and data analysis 

The method to extract the heavy-light decay constant is standard. We define a local and 
a smeared operator for the pseudoscalar channel by 

0^(t, x) = ^/(t, x)^5Mt, x), (13) 

oUt,x)=Y^Mt,xh,Mt,y)<i^''''''i\x-y\). (14) 

y 

in the Coulomb gauge. For the smearing function we use the form 0^^^(|a;|) = exp(— a|a;|''), 
with the parameters a and b chosen so as to reproduce the functional form of the heavy- 
light meson wave function measured in our simulations. We measure the two-point functions 
given by 

CWpitf,^) = Y.{O''pitf^^f)Op{U,0)), (15) 

Xf 

cfp{tf,u) = j:(oUtf,^f)ol\uM, (16) 

C%pitf,U) = T.iJilitf^^f)Ol\u,0)), (17) 

with the Dirichlet boundary condition in temporal direction. In this measurement the source 
is placed at the time slice tj=6 (at /?=5.7), 7 (5.9) and 16 (6.1). For the heavy-light meson 
with zero spatial momentum, C^pi)p(tj, tj) and C^p2)p(^/, ^i) are identical by construction. 
We fit the correlators to the exponential form, 

CWp{tf,U) ^ Zj^lexp{-aE'''-{tf - U)), (18) 
Cfp{tf,U) ^ Zfpexpi-aE'^'^tf - U)), (19) 
C%p{tf,td ^ Z%pexp{-aE'''-{tf-t,)), (20) 
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over a range of t where we find a plateau in the effective mass plot. Representative effective 
mass plots for Cpp(t/,tj), Cpp(t/,tj), Cjfo-,p(tf,ti), and Cjfi)p(tf,ti) are shown in Figures]^ 
(0) for the case of the heaviest (lightest) quark masses at (3=6.1. The signal is remarkably 
clean even for Cjfi^p which includes a spatial differential operator. To constrain the fit as 
much as possible we take the binding energy i?^" to be common among the correlators. This 
is particularly necessary for a stable extraction of Zpp since the signal for Cpp(t/, ti) is much 
noisier than those for the others. We estimate statistical errors of the fitted parameters using 
the jackknife method with unit bin size. Statistical correlation of data between different time 
slices or between different mass parameters is neglected in the fitting. 



C. Heavy- light Meson Mass 

We calculate the pseudoscalar meson mass aMp from a sum of the renormalized heavy 
quark mass and the binding energy E^^^ through the formula 

aMp = Z^aMo -E + aE^''', (21) 

where E is the energy shift and Zm the kinetic mass renormalization of the heavy quark. 

The one-loop calculation of E and Zm was first carried out by Davies and Thacker |^ 
and by Morningstar ^j. We repeat the calculation for our action NRQCD-I. We write the 
perturbative expansion of E, Z.^ and the wave function renormalization Z2h as 

E = a,A, (22) 

Zm = l + asB, (23) 

Z2H = l + a,C, (24) 

and list A, B, and C in Table |. 



D. Heavy-light Decay Constant 

The pseudoscalar meson decay constant fp is constructed from the contribution of each 
operator defined by 

a^'HU^.r ^ ^(0|i,'i\|P> ^ (25) 

where ^1 — 3K/4/tcrit represents the tadpole-improved wave function normalization factor 
for light quark. Including the one-loop corrections, the decay constant is given by 

a'/\fp^p) = [1 + aspfy/\fp^pY'^ +j:aj^a'/\fp^pf\ (26) 

We note that a^/'^{fp^ Mp)^^'> = a^/'^{fp^/ Mp)^"^^ holds in the rest frame of the heavy-light 
meson. 

In Figures | and |^ we show ai?*^™ and a^/^(/pA/Mp)W as a function oil/ k together with 
a linear (solid lines) and a quadratic (dotted lines) fit. We employ the linear fit for chiral 
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extrapolation since the difference between the hnear and quadratic fits are neghgibly small 
compared with the errors of the data. The linear fit is also used for an interpolation to the 
strange quark. The values of aE^^^ and a^^'^{fpy/MpY^^ at k = Kent as well as those at k^i 
and Ks2 extracted in this way are summarized in Tables |T|, |l^ and 0. 



One of the points we discuss in detail below is the effect of O(asaAQCD) improvement in 
the static limit. For this purpose we need to extract the decay constant in the static limit. 

Figure 1^ shows the dependence of (fpy/Mp)^^^ as a function of 1/Mp for each /? where Mp 
is calculated by the tree-level formula. The physical scale of lattice spacing a is determined 
from the string tension. We fit the mass dependence by 



\io) 




aMp {aMpf ^ 
We also fit the mass dependence of 2aMQa^/'^{fp^M^)^^'^ by 



(27) 



2aMoa^/^(/pv/Mp)« = 2aM,a"\fp^ Mpf^ ( 1 + + FIF^ ) • (28) 
^ ^ static V aMp [aMpY j 

Similarly, Figure % shows the 1/Mp dependence of the ratio -2Mo(/pa/Mp)^^V(/pv/^)^°^- 
The functional dependence can also be parameterized as 



(/pVMF)<"> (/pyMp)'°> 



static 




(29) 



The values of a^/2(/pVMF)^°^ I static and 2aMQa^/^{fp^M^)^^'^\^t^tic are given in Table and 
0, respectively. 

To obtain the B meson decay constant at the physical B meson mass, we fit the 1/Mp 
dependence of the renormalized quantity fp^/Mp in the renormalization group invariant 
form $p = [as(Mp)/as(Mp)]^/^^/pA/Mp instead of fitting the contribution of individual 
operators and summing the results. 



V. STATIC LIMIT 

We begin discussion of our results with the lattice spacing dependence of the decay con- 
stant in the static limit. This limit has the advantage that errors that depend on the heavy 
quark mass such as ©(osAqcd/^) vanish, and hence we can see the effect of O(asaAQCD) 
improvement more clearly. 

According to the discussion in Sec. the contribution of J^^^^ vanishes in the static 
limit. From Eq. (^, the matching relation in the static limit for the axial vector current is 
given by 



A, 



4cont 



1 + a 0^°^ ■ J"-"' ■ + a o^"''^a J^"'^^ (30) 

L T "sf'static "^static ' "sPstatic ""^static 5 K^'^J 



M) (disc) j-(disc) 



where Pstltic Jstltic the naive static limit (except anomalous dimension) of p^"* and 
>^iatl- p£ic and jJtS are defined as 
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(disc) 
r static 



a J, 



(disc) 



static 



lim p2V2aMo, 
lim 2aMoJiS. 

aAfo^oo 



(31) 
(32) 



The numerical value of the matching coefficients in the static limit is given in Table |. 
The decay constant is calculated from 



^static 



(33) 



with 



static 



II (0) 
1 + "sPstatic 



(/pjMp)(o) 



static 



(disc) 
'^sPstatic 



2aMo{fpJMp 



1(1) 



(34) 



A nominal value of Mo=4.5 GeV is used for the heavy quark mass to evaluate the logarithm of 
Pstatic- For the strong coupling constant we employ ay(g*) |]16| evolved from /i = 3.40/a 
to q*. Within one-loop calculations there is an uncertainty in the choice of the scale q*. 
We take the average of results obtained with q* = n/a and with 1/a, and consider the 
discrepancy from the two choices of q* as an upper and lower bound for the error due to 
two-loop corrections in the renormalization factor. 

Figure ^ shows the a dependence of the decay constant in the static limit f^'^^^'^ and 



j^tatic Qpgn symbols represent the results which are not corrected for the mixing effect of 



the operator a Jg^lu'c 



(2) 

which corresponds to the static limit of 2aMoJiJ^) , and filled symbols 
include this effect. Statistical errors are shown with solid bars, and uncertainties due to the 
choice of q* by dotted bars. From the figure we see that an apparent a dependence for the 
unimproved results is removed by the inclusion of the higher dimensional operator jft^ati'c 
the one-loop level. 

A worry with this observation is a sizable systematic error due to two-loop uncertainties. 
On this point we note that an estimate for the optimal value of q* for the multiplicative 
renormalization coefficient is known to be q* = 2.18/a for the combination of the static heavy 
quark and the unimproved Wilson light quark |2l|]. Since there seems to be no obvious reason 
that this value changes significantly for the 0(a)-improved light quark action, taking the 
difference of the results for q* = ir/a and 1/a may well be an overestimate of the two-loop 
uncertainty. An alternative estimate employing as{2/a) would reduce the error estimate by 
roughly a factor two. Furthermore the magnitude of this error is correlated among different 
/?, and between results without and with the improvement at each (3. Hence a reduction of 
the a dependence is less affected by the choice of q* than f^^^^'^ itself. 

Numerically the magnitude of the 0{asa) term relative to that of the leading operator 
Jg^itic (the static limit of j/att) takes values. 



(disc) . , . ifpVW)^''' 

X Pstatic X 2aMo 



(/PV^)(0) 

0.272(84) X 1.036 x (-0.712(9)) 
0.217(53) X 1.036 x (-0.622(7)) 
0.189(40) X 1.036 x (-0.546(8)) 



-0.201(84) at P=5.7, 
-0.140(53) at p=5.9, 
-0.107(41) at p=6.1, 



(35) 
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where the error is dominated by the uncertainty in as- At (3 = 6.1 the effect reduces f^^^^'^ 
by about 10% from the value without the improvement term. 

VI. B MESON DECAY CONSTANT 
A. Dependence on heavy-light meson mass 

In Figure | we present $p = [a^iMp) / as{MB)f'^^ fpy/M^ as a function of l/Mp for 
three values of (3 of our simulation. Open symbols denote results from the leading operator 
alone, and filled symbols show how they change due to the inclusion of the higher dimen- 
sional operators J^llt and J^Hf The factor [as{Mp) / UsiMB)]^^^^ is introduced to cancel the 
logarithmic divergence (l/vr) In(aMo) in the one-loop coefficient . For as{Mp) we use 
1^ evolved from /i = 3.40/a to Mp. The chiral limit is taken for the light quark. 
Solid and dotted error bars show the statistical error and the uncertainty due to two-loop 
corrections in the renormalization factors. The latter is estimated in the same way as for 
the static hmit discussed in Sec. 0. 

As first observed in Refs. [§,|l3,|ll|] , the contributions from the operators and 
sizably affects the decay constant. The dominant effect arises from J^J^^. A larger difference 
between the two sets of results toward the static limit is explained by the fact that the 
one- loop coefficient p^''/2aMo increases larger toward this limit (see Figure |l|). In contrast, 
the contribution of J^H^ is negligible since the perturbative coefficient p"^ /2aMQ stays very 
small (Ip^V^c-MoI <0.2) for our heavy quark mass aMo >1.2. 

As was the case for the decay constant in the static limit, uncertainties due to two- loop 
corrections are sizable, particularly at (3=5.7. This uncertainty does decrease, however, for 
weaker couplings at (3 = 5.9 and 6.1. It also becomes smaller as one moves down from the 
static limit toward the physical point for the B meson Mp = Mb- 

B. Dependence on lattice spacing 

By interpolating data shown in Figure ^ to the physical B meson mass, we obtain J'b for 
each p. The decay constant /b^ for Bg meson is calculated in a similar manner. The bare 
b quark mass that reproduces the physical B meson mass is listed in Table [V|, and and 
Jb^ at each (3 are given in Table [VI I| for the two choices of the scale q* = n/a and 1/a. 



The lattice spacing dependence of and /^^ is shown in Figure |Ty. Looking at the 
central values, we observe that a large a dependence exhibited by the results without the 
operator mixing (open symbols) is removed in the full result (filled symbols). This feature is 
clearer for /^j a variation is seen for J'b between j3 = 5.9 and 6.1, albeit with simultaneously 
larger statistical errors. While we need to keep in mind the uncertainty due to a choice of 
as, this result indicates that the lattice spacing dependence of the B meson decay constant 
is sizably reduced after including the 0{asa) and 0{as/M) mixing terms. 
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C. Estimate of systematic errors 



We now discuss possible sources of systematic errors and attempt an order estimate of 
their magnitude. 

As already discussed the uncertainty from the scale for the strong coupling constant, 
which is an 0{a1) effect, is quite significant. The magnitude of this error, estimated as half 
the difference of values for q* = n/a and 1/a is given in Table VIII for each /3. 

We employ a light quark action which is 0(a)-improved at one-loop level. Since the 
two-loop uncertainty in this improvement of ©(a^aAqco) is negligibly small, we expect the 
leading discretization error from the light quark sector to be ©((aAqcD)^), which is also the 
magnitude of scaling violation in the gluon sector. With a nominal value Aqcd= 300 MeV, 
we estimate its size to be 2%-8% depending on /3 as hsted in Table |V111 . 

Our results are obtained for NRQCD-I which represents the leading term in an expansion 
in 1/M. We examine possible corrections due to the truncation by explicitly comparing the 
results of NRQCD-I with those of NRQCD-II which is correct to 0{1/M'^). Figure [n| shows 
this comparison. We find that the correction does not exceed the statistical error, 

which is about 4% in the B meson mass region, as previously observed in Ref. [Q. Higher 
order uncertainties are expected to be even smaller. 

Another source of the systematic error is the perturbative matching of the action and 
the operators of NRQCD. In the one-loop calculation of the self-energy and the current 
renormalization, we have consistently included all terms of order 1/M. Hence 0{as/{aM)) 
corrections are properly taken into account in our calculation, and the leading error is of 
0(as/(aM)^). An order estimate for 0(as/(aM)^) is given in Table |VIII| for each (3. The 
magnitude increases for larger (3 since aM becomes smaller. 

Adding these four leading systematic errors in quadrature, we find the systematic error 
to be about 8% at [3=6.1 and 5.9, while it is much larger (~ 15%) at /5=5.7. We, therefore, 
take an average over the two results at P=6.1 and 5.9 to quote our final result. In doing so, 
we have confirmed that the difference between 6.1 and 5.9 is within the estimated systematic 
error, as it should be. 

In addition to the above systematic uncertainties, we must include an uncertainty in the 
lattice scale 1/a. Throughout this work we have used the scale set with the string tension 
^/o^. Taking a variation of the ratio rrip/ ^ over (5 = 5.9, 6.1 and 6.3, we assign a 3.5% error 
in the lattice scale as we did in Ref. [E3 



D. Results 

Our final result for the B meson decay constant in the quenched approximation is given 

by 

fs = 167(7)(15) MeV, (36) 
fs^ = 191(4)(17)(^^) MeV. (37) 

Here the central value is an average over the results at (3=6.1 and 5.9, and the errors are 
statistical and systematic in the given order. The systematic error includes 8% as estimated 
in the previous subsection and the error in the lattice scale of 3.5%, added in quadrature. 
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For Jb^ there is an additional uncertainty from the strange quark mass. We take the value 
from the K meson mass (ksi) for our central value. Employing the mass (^^2) gives a 
larger fs^, which is given in the final error. 



Our resuh is higher than that of Ali Khan et al. |TT| at /5=6.0 (/b = 147(11) (16) MeV and 
fBs= 175(8)(18) MeV), though consistent within one standard deviation. To compare with 
the results obtained with the the 0(a)-improved clover action, we quote from the studies of 
the Fermilab ||2^ group and of JLQCD 



= 164{tu){8) MeV (Fermilab), 
= 173(4)(13) MeV (JLQCD), 
/^^ = 185(;^^)(9) MeV (Fermilab), 
= 199(3)(14) MeV (JLQCD). 

Our values with NRQCD are in good agreement with these relativistic calculations. 



E. fsJfB 

Many of systematic uncertainties that appear in the calculation of the pseudoscalar decay 
constant /p(.^j cancel, if we consider the ratio fpjfp. In particular, the two-loop uncertainty 
in the matching of the axial current cancels out explicitly. 



Figure |T2| presents the 1/Mp dependence of fpJfp- We observe only a mild 1/Mp 
dependence. The difference between NRQCD-I and NRQCD-II is much smaller than the 
statistical error, which shows that the contribution of the 1/M^ terms is negligible. Finally, 
plotting the ratio as a function of lattice spacing (see Figure |l^), we find the results at three 
/3 values to be consistent with each other within the estimated errors. 

Our result is 

/B.//B = 1.15(3)(l)(^g), (38) 

which is obtained by averaging data at (3=6.1 and 5.9. The errors given are those from 
statistical, systematic and uncertainty in Kg. For the systematic error the leading contri- 
bution arises from the lack of the strange quark mass dependent renormalization, which 
is 0{asaAQCT)) for the deviation fss/ fs — 1- The magnitude is 3-5% assuming the order 
counting. 

VII. MASS SPLITTINGS 

A byproduct of our simulation is the mass difference between the B and Bg mesons, which 
can be compared with experiment. Since the heavy quark mass cancels in this difference, 
there are no direct perturbative corrections for this quantity, though they enter implicitly 
through determination of the bare 6-quark mass. 



We plot the 1/Mp dependence of the Bg — B mass difference in Figure |14|, where we 



observe the dependence to be small. The lattice spacing dependence is shown in Figure IT3. 



A variation of about 20%, beyond the statistical error of 8%, is seen between the data at 
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(3=6.1 and 5.9, which may represent scahng violation effects. Our result from an average 
over (3=6.1 and 5.9 is given by 



Mb, -Mb = 87(7) (4) (tf) MeV, (39) 

where, as before, the errors represent statistical, systematic and uncertainty. The possible 
systematic error is 0{{aAQCB)'^), which is 2-3%, and the uncertainty of l/a~ 3.5%. We 
estimate the error at 5% adding them in quadrature. The dominant error is the uncertainty 
of K,s- It is encouraging that our result agrees with experimental value 90 ± 2 MeV. 

The hyperfine splitting Mb* —Mb is another experimentally measured quantity. Previous 
lattice studies in the quenched approximation have found that the hyperfine splitting of 
heavy-light and heavy-heavy mesons are much smaller than experiment A possible 



reason for this discrepancy is an inappropriate value of the coupling ci for the gcr ■ B /2Mq 
term, for which we use the tadpole improved tree-level value. Since the hyperfine splitting of 
heavy-light mesons is proportional to Ci, and that of heavy-heavy mesons to c^, it is possible 
that large corrections of 0(0;^) remain (the non-perturbative calculation of this coupling 
has been done in Ref. which reports the possible 0{as) correction). Another possible 
source is the quenched approximation. 

The 1/Mp dependence of the hyperfine splitting obtained in our simulation is shown in 
Figure |16|. We observe that in the static limit the splitting linearly vanishes due to the heavy 



quark symmetry. Figure ^ shows the lattice spacing dependence of the splitting together 
with the experimental value of Mb* — Mb = 45.8 ± 0.4 MeV. While scaling is reasonably 
satisfied with our results, the magnitude is far below experiment. Averaging the two values 
at (3=6.1 and 5.9, we find 

Mb* -Mb = 23(7) (5) MeV, (40) 
Mb* - Mb, = 26(3)(5)(1J) MeV, (41) 

where we assume a 20% systematic error since the 0(0;^) correction is not known for Ci. 



VIII. CONCLUSIONS 

In this article we have presented a scaling study of the heavy-light meson decay constant 
using lattice NRQCD with the heavy-light current improved to 0{asa) and to 0{as/M) 
consistently to the one-loop order in perturbation theory taking into account mixings with 
the relevant higher dimensional operators. We have found the effect of the improvement to 
be substantial: the large a dependence of Jb is almost removed by the improvement. This 
is most apparent in the static limit, where the effect is purely 0{asa), but a similar effect is 
also seen for the physical B mass. 

The two main sources of systematic errors in our results are perturbative two-loop correc- 
tions, O(a^), in the renormalization factors for the NRQCD action, 0(a)-improved Wilson 
action and the axial vector current, and the one-loop corrections of 0(as/(aM)^) in the 
coefficients of the NRQCD action and the axial vector current. The former uncertainty is 
quite large at /5 = 5.7, but diminishes to a 5% level at weaker couplings of (3 = 5.9 and 6.1. 
The latter error, on the other hand, increases toward smaller lattice spacings, reaching ~ 
6% a.t (3 = 6.1. This counter increase of error represents the limitation of lattice NRQCD 
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itself; the method breaks down once the heavy quark mass becomes smaller than the inverse 
lattice spacing. 

The validity of a lattice NRQCD calculation of Jb hinges on the existence of a window 
in lattice spacing over which the two errors as well as scaling violations are small. We find 
these conditions to be satisfied at /3 = 5.9 — 6.1 where the combined systematic errors are 
estimated to be 8% (including the truncation in the 1/M expansion). Achieving better 
accuracy with NRQCD would require two-loop calculations to extend the window toward 
larger lattice spacings where the 0(as/(aM)^) error will become smaller. 

Another method for calculating heavy quark quantities on the lattice is the non- 
relativistic interpretation of relativistic actions [Q. The advantage is that a continuum 
extrapolation can be carried out. The simulations by Refs. p3| , p^^ have shown that the 
a dependence in the heavy-light decay constant is small for currently accessible range of 
f3 = 5.7 ~ 6.3. A continuum extrapolation, with either constant or linear fit in the lat- 
tice spacing, has yielded the decay constants with a systematic error of about 10%. A 
subtle point with this method, however, is that the a dependence of systematic errors is 
non-linear. Hence, strictly speaking, it is not correct to extrapolate the result with a simple 
linear or quadratic function of a. To achieve a prediction of the B meson decay constant 
more accurate than is available, one needs to improve the action and currents so that the 
systematic errors at finite values of /3 are further reduced. In this sense the study of O^aso) 
improvement should be necessary. 

In spite of the limitations still inherent in the present results with NRQCD and with 
relativistic actions, we find it encouraging that the efforts with the two approaches have 
now yielded predictions for the B meson decay constant in mutual agreement within the 
total error of about 10% for quenched QCD. 
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FIG. 1. l/aMo dependence of the one-loop coefficients for the axial vector current. Circles 
represent p^"* — l/7rln(oMo). Diamonds and triangles are p^'*/2aMo and /2aMQ respectively. 
The static limit is shown with the filled symbols. 
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FIG. 2. Effective mass of various correlators at /3=6.1 and (aMo, n) = (2.1, 2). The fitted value 
of aE^^^ is shown by a solid line, and the error is indicated by dashed lines. The light quark 
hopping parameter /c=0. 13586 is our heaviest one. 
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FIG. 3. Same as Figure 0, but with our lightest hght quark mass k=0. 13716. 
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FIG. 4. Chiral limit of the heavy-light binding energy aE^''' at (3=6.1 and (aMo,n) = (2.1,2). 
Open diamonds represent our data. Filled diamonds are the results in the chiral limit (Kcrit) or in 
the strange quark mass (ksi or ^^2) with linear fitting (solid line), and open squares are the results 
with quatratic fitting (dotted line). 
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FIG. 5. Chiral limit of the decay constant a3/2(/p^/Mp)(0) (left) and -2aMoa3/2(/p^/Mp)(i) 
(right) at /3=6.1 and (aMg, n) = (2.1, 2). The meaning of the symbols is the same as that in Fig. H, 
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FIG. 6. 1/Mp dependence of {fpy/ Mp)^^\ We used tree level value for Mp in the plot. Data 
at three /3 values are shown: /3=5.7 (diamonds), 5.9 (squares), and 6.1 (circles). The static limit 
(filled symbols) is obtained with a quadratic extrapolation. 
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FIG. 7. Ratio of the leading and mixing operators -2Mo(/p^/Mp)'^^V(/p\/Afp)^°^- We used 
tree level results for 1/Mp in the plot. Data at three (3 values are shown: /3=5.7 (diamonds), 
5.9 (squares), and 6.1 (circles). The static limit (filled symbols) is obtained with a quadratic 
extrapolation. 
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FIG. 8. The lattice spacing dependence of at k = /«crit (left) and k = Kgi (right). 

Open diamonds represent the result without the operator mixing, while filled diamonds include 
the mixing effect. The symbols show the q* averaged results, and are slightly shifted in horizontal 
axis so that error bars do not overlap. Solid error bars show the statistical error, and dashed ones 
show the uncertainty of q* from the difference of the two choices of q* = ir/a and 1/a. 
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FIG. 9. $p as a function of l/Mp for each /3. Open diamonds represent the result without 
the operator mixing, while filled diamonds include the mixing effect. The symbols show the q* 
averaged results. Solid error bars show the statistical error, and dashed ones show the uncertainty 
of q*. 
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FIG. 10. a dependence of fs (left) and fB^^iright). Filled symbols represent the result with the 
contribution from J^^^^ and j/^^t- Op^n symbols do not include these effects. The symbols show 
the q* averaged results, and are slightly shifted in horizontal axis so that error bars do not overlap. 
Solid error bars show the statistical error, and dashed ones show the uncertainty of q*. 
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FIG. 11. Comparison of $p from NRQCD-I (filled circles) with that from NRQCD-II (open 
circles) at tree level. 
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FIG. 12. 1/Mp dependence of fp^/fp with k^i. The symbols show the q* averaged results. 
Solid error bars show the statistical error, and dotted ones show the uncertainty of q*. 
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FIG. 13. a dependence of fpsl Ib with Ks\ (filled circles) and with Ks2 (open circles). The 
symbols show the q* averaged results. Solid error bars show the statistical error, and dotted ones 
show the uncertainty of q*. 
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FIG. 14. 1/Mp dependence of Mp^ — Mp with Kgi- We used the tree level results for 1/Mp in 
the plot. 
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FIG. 15. a dependence of Mb^ — Mb with (filled circles) and with Ks2 (open circles). The 
experimental value is shown by a solid line. 
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FIG. 16. 1/Mp dependence of Mp* — Mp. We used the tree level results for 1/Mp in the 
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TABLES 



aMo 


n 


(0) 

Pa - 


(l/vr) log(aMo) 


py /2aMo 


(2) /(^ jt r 

py /2aMo 


A 


B 


C 


oo 






-1.317 


0.000 


1.036 


1.069 




0.481 


12.0 


2 




-1.162 


0.026 


0.851 


1.022 


-0.091 


0.312 


10.0 


2 




-1.131 


0.030 


0.809 


1.011 


-0.075 


0.279 


7.0 


2 




-1.061 


0.036 


0.737 


0.983 


-0.032 


0.197 


6.5 


2 




-1.043 


0.037 


0.725 


0.976 


-0.023 


0.177 


5.0 


2 




-0.970 


0.040 


0.656 


0.946 


0.018 


0.094 


4.5 


2 




-0.937 


0.039 


0.628 


0.931 


0.036 


0.055 


3.8 


2 




-0.876 


0.036 


0.578 


0.903 


0.069 


-0.014 


3.5 


2 




-0.846 


0.034 


0.559 


0.888 


0.086 


-0.052 


3.0 


2 




-0.782 


0.026 


0.516 


0.855 


0.119 


-0.127 


2.1 


3 




-0.626 


-0.015 


0.442 


0.754 


0.329 


-0.315 


1.5 


3 




-0.433 


-0.108 


0.378 


0.621 


0.456 


-0.542 


1.3 


3 




-0.341 


-0.173 


0.360 


0.547 


0.511 


-0.647 


0.9 


4 




-0.088 


-0.445 


0.374 


0.300 


0.782 


-0.921 



TABLE L One- loop coefficients of the axial vector current p^\ p^ and p)^ defined in Eq. @. 
The self-energy corrections are also listed. 





6.1 


5.9 


5.7 


size 


24^ X 64 


16^ X 48 


12^ X 32 


T^conf 


120 


300 


300 




1.525 


1.580 


1.674 


K 


0.13586 


0.13630 


0.13690 




0.13642 


0.13711 


0.13760 




0.13684 


0.13769 


0.13840 




0.13716 


0.13816 


0.13920 


Uq 


0.8816 


0.8734 


0.86087 


(aMo,n) 


(7.0,2) 


(10.0,2) 


(12.0,2) 




(3.5,2) 


(5.0,2) 


(6.5,2) 




(2.1,2) 


(3.0,2) 


(4.5,2) 




(1.5,3) 


(2.1,3) 


(3.8,2) 




(0.9,4) 


(1.3,3) 


(3.0,2) 




0.149 


0.164 


0.188 


ay(l/a) 


0.229 


0.270 


0.355 




0.13767 


0.13901 


0.14157 


Ksl 


0.13635 


0.13702 


0.13800 




0.13609 


0.13657 


0.13707 


1/a (GeV) 


2.29 


1.64 


1.08 



TABLE II. Lattice parameters. 
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aMo 




q* = -K/a 


= 1/a 




q* = vr/a 


q* = 1/a 


aE^'"^ 


q* = n/a 


q* = 1/a 












=5.7 










12.0 


0.669(11) 


12.271(11) 


11.919(11) 


0.737(7) 


12.339(7) 


11.987(7) 


0.755(7) 


12.357(7) 


12.005(7) 


6.5 


0.670(10) 


6.958(10) 


6.770(10) 


0.739(6) 


7.027(6) 


6.839(6) 


0.758(5) 


7.046(5) 


6.858(5) 


4.5 


0.665(8) 


5.021(8) 


4.892(8) 


0.738(5) 


5.093(5) 


4.965(5) 


0.758(4) 


5.113(4) 


4.985(4) 


3.8 


0.663(8) 


4.341(8) 


4.235(8) 


0.737(5) 


4.416(5) 


4.309(5) 


0.757(4) 


4.436(4) 


4.330(4) 


3.0 


0.658(7) 


3.564(7) 


3.481(7) 


0.735(4) 


3.642(4) 


3.559(4) 


0.756(4) 


3.663(4) 


3.580(4) 












=5.9 










10.0 


0.531(8) 


10.244(8) 


10.058(8) 


0.580(5) 


10.292(5) 


10.106(5) 


0.591(5) 


10.303(5) 


10.117(5) 


5.0 


0.528(7) 


5.389(7) 


5.298(7) 


0.575(4) 


5.435(4) 


5.345(4) 


0.586(4) 


5.446(4) 


5.356(4) 


3.0 


0.522(6) 


3.440(6) 


3.387(6) 


0.569(3) 


3.488(3) 


3.435(3) 


0.580(3) 


3.498(3) 


3.446(3) 


2.1 


0.511(6) 


2.601(6) 


2.594(6) 


0.560(3) 


2.650(3) 


2.643(3) 


0.571(3) 


2.661(3) 


2.654(3) 


1.3 


0.487(5) 


1.806(5) 


1.818(5) 


0.539(3) 


1.859(3) 


1.871(3) 


0.552(2) 


1.871(2) 


1.883(2) 












=G.l 










7.0 


0.435(8) 


7.255(8) 


7.158(8) 


0.471(5) 


7.292(5) 


7.195(5) 


0.479(5) 


7.299(5) 


7.202(5) 


3.5 


0.423(7) 


3.835(7) 


3.788(7) 


0.462(4) 


3.875(4) 


3.828(4) 


0.470(4) 


3.883(4) 


3.836(4) 


2.1 


0.408(6) 


2.499(6) 


2.494(6) 


0.450(3) 


2.541(3) 


2.535(3) 


0.458(3) 


2.549(3) 


2.544(3) 


1.5 


0.393(5) 


1.903(5) 


1.908(5) 


0.436(3) 


1.946(3) 


1.951(3) 


0.445(3) 


1.954(3) 


1.959(3) 


0.9 


0.359(4) 


1.320(4) 


1.352(4) 


0.404(3) 


1.364(3) 


1.396(3) 


0.413(2) 


1.373(3) 


1.405(2) 



TABLE III. Binding energy and the total mass of the heavy-hght pseudoscalar mesons. 



28 



aiviQ 


K — Hcvit 


1^ — l^sl 


K = Ks2 








5.7 






oo 


U. 575(41 j 






0.814(34) 


0.851(36) 


ion 


0.5oo(zo j 






0.693(19) 


n '700/on\ 
U.7zz(zUj 


D.O 


U.OoHiy j 






0.615(13) 


U.Doo(io j 


4.0 


U.4oi(i0 j 






0.556(11) 


U.O /0(iU j 


o.o 


U. 400(^14 J 






0.527(9) 


u.o4D(y j 


o.U 


U.4zl(^lz J 






0.486(8) 


U.OUo(o J 






13 = 


5.9 






OO 


n oi 0/1 c;^ 
O.olz(15j 






0.370(11) 


r\ OOO /1 1 \ 

0.383(11 j 


inn 
iU.U 


n OQc^/' 1 1 ^ 
U.zoO^ii j 






0.333(8) 


U.o44( r j 


o.U 


U.zDU(y j 






0.296(7) 


U.oU4( / j 


o.U 


U.zoo^o j 






0.264(5) 


n OV1 
U.Z I i(0 J 


1 








0.240(4) 


U.Z4D(4 ) 


1 Q 

i.o 


U.i (c>\p ) 






0.201(3) 


U.ZU (\o) 






P = 


^ 1 

D.l 






OO 


0.178(12) 






0.205(9) 


0.210(8) 


7.0 


0.159(9) 






0.185(7) 


0.190(6) 


3.5 


0.140(7) 






0.165(5) 


0.170(4) 


2.1 


0.124(5) 






0.148(4) 


0.152(3) 


1.5 


0.114(4) 






0.135(3) 


0.140(3) 


0.9 


0.096(3) 






0.114(2) 


0.118(2) 



TABLE IV. Raw data a^l'^{fp^fM^)^^^ at Kcrit, «si and Ks2- 
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aMo 




1^ — l^sl 


K = Ks2 


P=5.7 


oo 


-0.485(34) 


-0.556(27) 


-0.576(29) 


12.0 


-0.455(22) 


-0.511(16) 


-0.526(16) 


6.5 


-0.436(18) 


-0.482(12) 


-0.495(12) 


4.5 


-0.415(15) 


-0.458(10) 


-0.470(10) 


3.8 


-0.403(14) 


-0.446(9) 


-0.458(9) 


3.0 


-0.387(12) 


-0.429(8) 


-0.441(8) 



/3 = 5.9 



oo 


-0.194(12) 


-0.226(8) 


-0.234(8) 


10.0 


-0.189(8) 


-0.215(6) 


-0.221(5) 


5.0 


-0.183(7) 


-0.203(5) 


-0.208(5) 


3.0 


-0.176(7) 


-0.193(4) 


-0.198(4) 


2.1 


-0.170(6) 


-0.187(4) 


-0.191(4) 


1.3 


-0.158(5) 


-0.176(3) 


-0.180(3) 


^ = 6.1 


oo 


-0.098(8) 


-0.111(6) 


-0.113(5) 


7.0 


-0.092(6) 


-0.105(4) 


-0.108(4) 


3.5 


-0.086(5) 


-0.100(3) 


-0.103(3) 


2.1 


-0.082(4) 


-0.097(3) 


-0.100(3) 


1.5 


-0.080(4) 


-0.095(2) 


-0.098(2) 


0.9 


-0.079(3) 


-0.093(2) 


-0.096(2) 



TABLE V. Raw data of 2aMoa^/2(/pVMp)(^) at Kcrit, ^si and Ks2- 
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l3 = Q.l 


P = 5.9 


/3 = 5.7 


tree 


1.903(5) 


2.710(6) 


4.206(8) 


q* = n/a 


1.913(6) 


2.765(6) 


4.345(9) 


q* = 1/a 


1.919(6) 


2.804(7) 


4.476(9) 



TABLE VI. Bare b-quark mass that reproduces the physical B meson mass. 





/3 = 6.1 


13 = 5.9 


/3 = 5.7 




tree (f = Tx/a q* = 1/a 


tree g* = vr/a q* = l/a 


tree q* = ii/a q* = 1/a 


Ib 

fBs {l^sl) 
fBs {l^s2) 


0.184(7) 0.166(7) 0.157(6) 
0.215(5) 0.195(5) 0.185(4) 
0.222(5) 0.201(4) 0.190(4) 


0.210(6) 0.181(6) 0.163(5) 
0.233(4) 0.203(4) 0.183(3) 
0.239(4) 0.208(4) 0.187(3) 


0.233(7) 0.189(6) 0.150(5) 
0.265(5) 0.217(4) 0.174(3) 
0.274(5) 0.225(4) 0.180(3) 



TABLE VII. Results for fe and /b, in GeV. 





13 = 6.1 


P = 5.9 


(3 = 5.7 




3% 


5% 


12% 


OiiaAQCD?) 


2% 


3% 


8% 




6% 


•1% 


2% 



TABLE VIII. An order estimate of the possible systematic errors. 
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